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CZ2 

C3 ABSTRACT 



Context. The existence of debris disks around old main sequence stars is usually explained by continuous replenishment of small dust 
grains through collisions from a reservoir of larger objects. 

Aims. We present photometric data of debris disks around HIP 103389 (HD 199260), HIP 107350 (HNPeg, HD 206860), and 

HIP 1 14948 (HD 219482), obtained in the context of our Herschel Open Time Key Program DUNES (DUst around NEarby Stars). 

Methods. We used Herschel/PACS to detect the thermal emission of the three debris disks with a 3cr sensitivity of a few mJy at 

100 //m and 160 //m. In addition, we obtained Her.schel/PACS photometric data at 70//m for HIP 103389. These observations are 

complemented by a large variety of optical to far-infrared photometric data. Two different approaches are applied to reduce the 

C^ Herschel data to investigate the impact of data reduction on the photometry. We fit analytical models to the available spectral energy 

f^ distribution (SED) data using the fitting method of simulated thermal annealing as well as a classical grid search method. 

04 Results. The SEDs of the three disks potentially exhibit an unusually steep decrease at wavelengths > 70/jm. We investigate the 

T-H significance of the peculiar shape of these SEDs and the impact on models of the disks provided it is real. Using grain compositions 

IL" that have been applied successfully for modeling of many other debris disks, our modeling reveals that such a steep decrease of the 

. 5^ SEDs in the long wavelength regime is inconsistent with a power-law exponent of the grain size distribution -3.5 expected from a 

\^ standard equilibrium collisional cascade. In contrast, a steep grain size distribution or, alternatively an upper grain size in the range 

%^ of few tens of micrometers are implied. This suggests that a very distinct range of grain sizes would dominate the thermal emission 

^ of such disks. However, we demonstrate that the understanding of the data of faint sources obtained with Herschel is still incomplete 

and that the significance of our results depends on the version of the data reduction pipeline used. 

Conclusions. A new mechanism to produce the dust in the presented debris disks, deviations from the conditions required for a stan- 
dard equilibrium collisional cascade (grain size exponent of -3.5), and/or significantly different dust properties would be necessary 
to explain the potentially steep SED shape of the three debris disks presented. 

Key words. Stars: circumstellar matter - Stars: individual: HIP 103389, HIP 107350, HIP 114948 - Infrared: planetary systems - 
Infrared: stars 

1. Introduction 

* Herschel is an ESA space observatory with science instruments 
provided by European-led Principal Investigator consortia and with im- Debris disks were first discovered via infrared excess emission 
portant participation from NASA. associated with main sequence stars such as Vega detected by 
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the Infrared Astronomical Satellite (IRAS; Aumann et al.|1984) l. 
The first spatially resolved image of a debris disk was that of 
/?Pictoris ( [Smith & Terrilej|1984 ) in optical scattered light. In 
the last few years, the Spitzer Space Telescope has revealed 
that debris disks are common around main sequence stars (e.g., 
Trilhng et al. 2008|l. Several Key Programs on the Herschel 



Space Observatory (Pilbratt et al. 2010 1 are dedicated to the 
study of various aspects of the formation and evolution of plane- 
tary systems and their attendant circumstellar debris disks (e.g., 
[Augereau et al. 2008 1. Our Herschel Open Time Key Program 
(OTKP) DUNES (DUst aiound NEarby Stars; |Eiroaet al.|2010| l 
aims to detect debris disks with fractional luminosities simi- 
lar to the Edgeworth-Kuiper B elt level (Ld/i-* = 10"^ to 10"''; 
Stern|1996 Vitense et al. 2010[) around a volume limited sample 



{d < 20 pc) of Sun-like stars (F, G, and K spectral type). Some 
additional sources at 20 pc < li < 25 pc are included because of 
their known excesses that have previously been detected with 
Spitzer, or because they are known exoplanet host stars. 

As the most readily detectable signposts of other planetary 
systems, debris disks help us to improve our understanding of 
the formation and evolution of them as well as of our own solar 
system ( [Meyer et al.|2007[|^^t|2008[[Kriv"ov|2010| l. Studying 
the spectral energy distribution (SED) of the dust alone usually 
provides only weak, ambiguous constraints to their properties 
such as chemical composition, grain size, and spatial distribu- 
tion ( [Wolf & Hillenbrand[[2003| l. For example, the location of 
the inner disk radius of the dust distribution is strongly degener- 
ate with the lower limit of the grain size distribution. However, 
the grain size distribution can be described by a power-law and 
the power-law exponent derived from pure SED fitting is usu- 
ally well constrained and consistent with the analytical value of 
-3.5 derived by Dohnanyi] (T969) under several assumptions: (1) 
The particles are produced through a coUisional cascade from 
infinitely large to extremely small grains. (2) The grains are not 
affected by any other influence such as stellar radiation and drag 
forces. (3) The strength of the particles, i.e., the energy per vol- 
ume fraction necessary to disrupt them, is independent from the 
grain size. This is further referred to as a standard equilibrium 
collisional cascade. 

In this paper, the potential //ersc/zeZ/DUNES discov- 
ery of an unusually steep decrease of the SEDs of three 
spatially unresolved debris disks around Sun-like stars, 
HIP 103389 (HD 199260), HIP 107350 (HNPeg, HD 206860), 
and HIP 114948 (HD 219482), is presentecQ It is demonstrated 
on the example of the three targets that the analysis of data of 
faint point sources obtained with Herschel depends very much 
on the photometric calibration and the exact determination of 
the uncertainties. In particular, the significance of our results 
depends on the version of the data reduction pipeline used. 
Provided it is real, the steep decrease occurs in the range of 
70/im to 160 ;um, inconsistent with a dust grain size distribution 
following a power-law derived from a standard equilibrium col- 
lisional cascade. The 100 ;um and 160 ;um data represent the first 
photometric measurements of these faint disks at wavelengths 
> 70/im and, thus, are the first ones to reveal the unusual shape 
of the SEDs. Results from detailed SED modeling including ad- 
ditional photometric data from the literature are presented as 
well. 



' Two of the sources, HIP 103389 and HIP 107350, are shared tai'gets 
between the DUNES survey in the context of which the analysis is car- 
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Fig. 1. Images of the three sources at 100 //m (top) and 160 //m 
{bottom) from Reduction 1 (Sect. |2|i. The images are displayed 
in a linear stretch from -O.lFpeak to Fpeak of the source of 
interest. The peak flux is 0.38 mJy /pixel, 0.23 mJy /pixel, and 
0.54 mJy/pixel at 100 /vm and 0. 17 mJy /pixel, 0.16 mJy /pixel, 
and 0.45 mJy /pixel at 160 yum for HIP 103389, HIP 107350, and 
HIP 114948, respectively. The white circles in the 160yL(m im- 
ages are centered on the brightest pixel of the source, respec- 
tively. 

Observations and data reduction are described in Sect.[2|to- 
gether with the basic observational results. A theoretical discus- 
sion on how steep the decrease of the SED of a debris disk to- 
ward longer wavelengths is expected to be and a characteriza- 
tion of unusually steep SED sources are given in Sect l3] A de- 
tailed description of our modeling of the systems can be found 
in Sect. [4| and results are discussed in Sect. [5] Conclusions are 
given in Sect.l6] 

2. Observations & data reduction 

Two Herschel/PACS mini-scan map observations of each target 
were taken with the 100/160 channel combination at array ori- 
entation angles of 70° and 110° providing scan and cross-scan 
coverage to assist in the removal of noise artifacts from the fi- 
nal composite mosaic. In addition, two scan map observations 
of HIP 103389 were taken with the 70/160 channel combina- 
tion with the same array orientations. Each scan map consists 
of 10 legs of 3' length, with a 4" separation between legs, at the 
medium slew speed (20" per second). In this way, a region of 
» 1 square arc minute around the source position was covered 
to uniform depth in the resulting mosaic. A summary of the ob- 
serving log is presented in Table [T] 

PACS data reduction was carried out in version 4.2 
{Reduction 7) and 7.2 {Reduction 2) of HIRE ( [Ott et al.[[20T0l l 
starting from the level products using a modified version of 
the standard reduction script (provided within HIPE). Two dif- 
ferent reductions have been carried out, because changes in the 
still developing pipeline of HIPE may significantly affect the re- 
sults. The separate scans at the two position angles of each chan- 
nel pair were mosaicked to produce a final image at each wave- 
length. Output scales for the final mosaics are 1 " per pixel for the 
70/jmand 100 ;um images and 2" per pixel for the 160 jum image 
(smaller than the native detector pixel sizes of 3'.'2 at 70 //m and 
100 jum and 6'.'4 at I60fj.m, respectively). 

For the further analysis of the data, different approaches have 
been used for Reductions 1 and 2. These are described in the 
following: 
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Table 1. Summary of HerschelfPACS observations of the three 
objects 



Object 


Observation ID 


OD 


A[fim] 


T[s] 


HIP 103389 


1342193157/58 


322 


70/160 


180 


HIP 103389 


1342193159/60 


322 


100/160 


1440 


HIP 107350 


1342195779/80 


355 


100/160 


1440 


HIP 114948 


1342196803/04 


372 


100/160 


1440 



Notes. OD is the observing date that counts the elapsed time in days 
since launch on 14 May 2009, A gives the central wavelength of the 
two filters used simultaneously, and T is the total on-source integration 
time. 



Reduction 1: A high-pass fiher was used to remove large scale 
background variation from the images, with filter widths of 3 1 
frames at 70;um and 100//m and 51 frames at 160//m (equiva- 
lent to ~ 60" and 100", respectively). A central region of 30" 
radius in the images was masked from the high pass filter pro- 
cess to prevent the removal of any faint, extended structure near 
the source position. PACS fluxes and sky noise were measured 
by aperture photometry using a custom script based on the IDL 
APER routine, based on the DAOPHOT aperture photometry 
module and cross-checked with the internal HIPE aperture pho- 
tometry routines. The aperture radius and sky annulus dimen- 
sions were 20" and 30" - 40", respectively at all bands. Sky 
noise for each wavelength was calculated from the rms vari- 
ance of the sky annulus pixel values multiplied by the beam 
size. Results were scaled by the appropriate aperture correction. 
Absolute uncertainties are 5% at 70 /im and 100 /im and 10% at 
160yum using the calibration files provided along with this HIPE 
version. 

The images of the three sources from Reduction 1 at lOOyum 
and 160 //m are shown in Fig.fT] 

Reduction 2: The same filter width as in Reduction 1 has been 
used for the high-pass filtering, but sources were masked based 
on a threshold criterion, the value of which was calculated from 
the standard deviation of all non-zero pixels in the image. PACS 
fluxes and sky noise were measured by aperture photometry us- 
ing both the internal HIPE aperture photometry routines and 
the MIDAS data analysis package to check for consistency. The 
source flux was measured within a circular aperture of 4", 5", 
and 8" radius at 70 /im, lOOyum, and 160/vm, respectively, and 
scaled by the appropriate aperture correction. Sky noise for each 
wavelength was calculated from the rms variance often sky aper- 
tures of the same size as the source aperture and randomly dis- 
tributed across the uniformly covered part of the image. The re- 
sulting deviation of these apertures was calibrated by the aper- 
ture correction and correlated noise correction factors. Extensive 
testing has shown that this approach efficiently accounts for the 
correlated noise in the output images. In fact, this approach pro- 
vides slightly conservative uncertainties. An extended discus- 
sion will be included in the survey overview paper (Eiroa et al., 
in prep.). Aperture correction factors have also been applied to 
the rms aperture flux dispersion. Photometric calibration uncer- 
tainties in this version of HIPE are 3% - 59^ 



of HIP 103389, HIP 107350, and HIP 114948 (Eiroa et al., in 
prep.). The stellar contribution to the total SED of each star 
is estimated using a PHOENIX/GAIA synthetic stellar model 
(Brott & Hauschildt 2005). Specific models for T^ff, logg* and 
[Fe/H] of each star are built by interpolation from the available 
grid. A normalization is applied to the short wavelength range 
{A < 20yum) of the long wavelength section of the Spitzer/]RS 
spectrum where available. Otherwise, the normalization is done 
to the flux measurements at photometric bands from B to K. 
All three sources are found to have significant excess (> 3 cr) 
at 70 /im and 100 //m. HIP 114948 also has significant excess at 
160 jum. HIP 103389 and HIP 107350 have estimated excesses 
at 160/im between 0.8cr and 2.4cr, depending on the reduc- 
tion. Offsets between optical positions and positions in the PACS 
70;um images are below 1 cr of the Herschel pointing accuracy. 
In each case, the observed excess is attributed to the presence of 
a debris disk associated with the star. The measured FWHM of 
the sources (derived from a 2-D Gaussian fit) are consistent with 
unresolved objects. One can therefore constrain the extension of 
the emitting area (adopting the distance of the objects listed in 
Table[3]l to a diameter of less than 7'.' 1 (156 AU) for HIP 103389, 
7'/l (126 AU)forHIP107350, and7'.'2(148AU) for HIP 114948. 
The spectral slopes of the three systems form Reduction 1 
and Reduction 2 are listed in Table |2l The SEDs of all three 
sources are found to exhibit an unusually steep decrease at the 
wavelength range of 70jum to 160/vm (steeper than a black 
body radiator in the Rayleigh- Jeans regime) from the results of 
Reduction 1 . The steepness is evaluated based on the spectral in- 
dex of the SEDs between two different wavelengths (see Sect. l3] 
for a detailed description). The significance of this steepness 
is evaluated using error propagation. However, the results from 
the new (but not necessarily better) version of the data reduc- 
tion (Reduction 2) give a shallower decrease of the SED toward 
longer wavelengths. Thus, we emphasize that this is a potential 
discovery of a steep decrease of the SEDs of the three sources 
and that a deeper understanding of the data (in particular for very 
faint sources) than possible until now is necessary to make a fi- 
nal statement. In the subsequent sections, we discuss the conse- 
quences of this steepness of the SEDs under the assumption that 
it is real. 



3. What is an unusually steep SED? 

Infrared excesses observed from debris disk systems are the re- 
sult of absorption of incident starlight by dust grains in the disk 
and subsequent re-emission of stellar radiation. Since the dust 
is colder than the star, its emission peaks at longer wavelengths 
(we assume that there is only one peak). In this section, a quali- 
tative discussion is given of how steep the SED of a debris disk is 
expected to be. We first investigate the shape of the excess treat- 
ing both the star and the disk as single temperature black body 
radiators. Later, we will discuss the effects of a more realistic 
disk model, i.e., allowing for a range of temperatures and more 
realistic grain properties. 



2.1. Observational results 

The data obtained in this work for the three sources using 
both reductions as well as the predicted photospheric fluxes 
are listed in Table l2] Table l3] lists the main stellar parameters 



Release note PICC-ME-TN-037 



3.1. Treating ttie dust as a single temperature black body 
The spectral index A of an SED is defined as 

d log Fy 



51ogi 



(1) 
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Table 2. Observational results obtained in this work" 



Source 



HIP 103389 



HIP 107350 



HIP 1 14948 



Reduction 1 Reduction 2 Reduction 1 Reduction 2 Reduction 1 



Reduction 2 



FpAcs70 [mjy] 


47.4 ±2.7 44.0 ±2.3 






^MiPS70 [mJy] 


46.6 ±3.8 


28.4 ±2.5 


68.7 ±3.0 


f PAcs 100 [mJy] 


23.7 ±1.4 26.3 ±1.7 


11.0±0.9 15.1 ±1.3 


42.5 ±2.2 40.8 ±1.6 


FpAcs 160 [mJy] 


5.0 ±1.3 7.7 ±2.5 


4.4 ±1.5 4.4 ±2.3 


12.7 ±1.9 13.3 ±2.2 


F*.i(i^m" [mJy] 


13.8 


13.4 


15.0 


f *,ioo^ra* [miy] 


6.8 


6.6 


7.3 


^*,i60,im* [mIy] 


2.7 


2.6 


2.9 


Al00,70 


1.94 ±0.32 1.44 ±0.33 


2.66 ±0.45 1.77 ±0.49 


1.35 ±0.26 1.46 ±0.23 


Al60,100 


3.31 ±0.69 2.61 ±0.82 


1.95 ±0.88 2.62 ±1.30 


2.57 ± 0.44 2.38 ± 0.43 


A 160,70 


2.72 ±0.39 2.11 ±0.46 


2.25 ±0.53 2.26 ±0.74 


2.04 ±0.22 1.99 ±0.25 



FWHMpAcsioo^ 

RApACS 100 

DECpAcs 100 

^^optical 



DEC 



optical 



6'.'8x7'.'l 
20"56™47:47 
-26°17'45'.'5 
20''56™47:33 
-26°17'47'.'0 



7'.'lx6'.'0 
2lH4"31!27 
+ 14°46'19'.'4 
2l''44"31!33 
+ 14°46'19'.'0 



7'.' 2 X 7'.' 2 
23"16"57»46 
-62°00'05'.'7 
23''16"57^:69 
-62°00'04'.'3 



Notes. Uncertainties of the flux measurements are total uncertainties including sky noise and calibration uncertainties. A^j „ is the spectral index 
measured from vi and V2 (SectlSl, identified with the corresponding PACS bands. 

*°' The SpitzerfMlPS 70 /im photometry is also listed, since these data complement the Herschel observations and are important to illustrate 
the unusual behavior of the SEDs. These data are pubhshed by |Beichman et al.|p006l HIP 103389 and HIP 114948) and Bryden et al. (2006] 
HIP 107350) and have been re-reduced in the context of DUNES (Eiroa et al., in prep.). ^ ' Predicted stellar photosphere using a PHOENIX/GAlA 
synthetic stellar model (Eiroa et al., in prep.). '"^^ FWHM as measured from the PACS images at 100 yum using a 2-D Gaussian fit. 



where Fy is the total flux of the object at the frequency v. It is 
typically measured as the slope of the SED between two data 
points available, i.e.. 



log Fy, - log Fy 

log V2 - log Vi 



(2) 



We further denote the two frequencies v\ and V2 with their cor- 
responding PACS wave bands (e.g., Aioo.vo)- The spectral index 
of a black body radiator amounts to A = 2 in the Rayleigh- Jeans 
regime. It is decreasing toward shorter wavelengths where the 
Rayleigh- Jeans approximation is not valid. A spectral index of 
A > 2 means that the SED falls off toward longer wavelengths 
steeper than a single temperature black body in its Rayleigh- 
Jeans regime. In the wavelength range where dust reemission 
is observed, the emission of the star can be very well approxi- 
mated by a black body radiator in its Rayleigh- Jeans regime. It 
will only contribute as an additional component with a spectral 
index of A = 2. Thus, the contribution of the star will - if sig- 
nificant - move the spectral slope of the whole system closer to 
2. 

3.2. More exact treatment of the disk 



Our first step to expand the above discussion is to allow for dif- 
ferent dust temperatures, while the dust is still treated as a num- 
ber of single temperature black bodies (i.e., a radially extended 
disk). A two component black body has been used successfully 
to model SEDs of debris disks (e.g., Hillenbrand et al. 2008 1. 
This results in a flatter SED (Av,_vj < 2) than in the above case, 
since the colder dust component will contribute additional flux 
in the long wavelength regime of the warmer dust component. 
Thus, using several components with different temperatures one 
can only reach Av.,,v, = 2 in the long wavelength regime of the 
coldest dust component (i.e., at even longer wavelengths com- 
pared to the single black body case above). At shorter wave- 
lengths Av,,v, will always be < 2. 
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Fig. 2. Absorption efliciency QibsW for astronomical silicate 



grains of different grain radii a (see Sect. 3.2 for details) 



Using Mie theory, one can derive absorption efficiencies 
Gabs (A) for more realistic dust grains. The emission of a dust 
grain is then described by the following equation which is a mod- 
ification of Planck's law: 



B'a (^, Tdust) 



1 



2hc^ 



A^ exp [hc/AkTaust] - 1 



gabs U) 



(3) 



The shape of the simple black body SED is modified by the 
function Qabs (^) which will result in a steeper decrease with 
wavelength in the Rayleigh- Jeans regime of the black body if 
2abs W is decreasing with wavelength. Fig.l2]shows gabs i^) for 
astronomical silicate ( Draine|[2003| l and different grain radii a. 
6abs (A) is close to 1 at short wavelengths, while it indeed ex- 
hibits a break atAx 2na and then decreases toward longer wave- 
lengths with A^^ (or v^). This is analog to the discussion of the 
opacity index, e.g., by |Draine| ( |2006[ ). An SED of a dust disk that 
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Table 3. Physical properties of the three stars considered in this work (for details and references see Sect. 



4.1: 



Object 



HIP 103389 



HIP 107350 



HIP 114948 



Name 



HR8013,HD 199260 



HN Pes, HD 206860 



HR 8843, HD 219482 



Distance (pc) 


21.97 


17.88 




20.54 


Spectral type and luminosity class 


F7V 


GOV 




F7 V 


Range of spectral types 


F6V-F8V 


GOIV- 


-V 


F6V-F8V 


V,B-V 


5.70,0.51 


5.96, 0.59 


5.64, 0.52 


Absolute magnitude My, bolometric correction 


3.99, -0.02 


4.70, - 


-0.05 


4.08, -0.02 


Bolometric luminosity, L^ [Lq] 


2.025 


1.090 




1.867 


Effective temperature (K) 


6257 


5952 




6240 


Surface gravity, log gt. 


4.36 


4.44 




4.31 


Radius, R^ [Ra\ 


1.2 


0.99 




1.17 


Metallicity [Fe/H] 


-0.14 


-0.07 




-0.21 


Radial velocity (km/s) 


-16.1 


18.1 




0.47 


Rotational velocity, v sin / (km/s) 


13.7 


12.8 




9.0 


Rotation period, P (days) 


4.4 (P/ sin 


4.7 




- 


Li I equivalent width (mA) 


80.0 


91.9 




73.7 


Space velocities U, V, W (km/s) 


-17.2,-10.2, 1.1 


-13.9, 


-20.1,-11.1 


14.9, -8.0, -5.7 


Activity, log i?[j^ 


-4.402 


-4.48 




-4.434 


X-Ray luminosity, logLx/L* 


-4.71 


-4.37 




-4.43 


Mass, M* [Mo] 


1.28, 


0.98 




1.02 


Age (Gyr) 


0.304, 0.184, 0.748(sin0 


0.530, 


0.162,0.291 


0.385,0.115,- 



Notes. The ages given are derived using the following methods (left to right): logR[^^ activity index levels, ROSAT X-Ray luminosities Lx/Lboi, 
Rotation period. For more details and references, see Sect.|4.1| 



is steeper in the long wavelength regime than a black body is ob- 
viously possible. However, debris disks are usually expected to 
be radially extended (e.g., Backman et al. 2009 Ertel et al. 2011 
Lohne et al.|2011[l and to contain particles of different size. This 



suggests a broad range of dust temperatures to be present result- 
ing in a flattening of the SED. In addition, the grain size depen- 
dence of the break in Qabs (/I) means that the grains have to be 
smaller than »; A/2n to result in an SED that falls off steeper than 
a black body ( |Draine|[2006) ). This suggests grains smaller than 
~ 15/im for the three disks presented in this paper (Table |2]i. 
Thus, one qualitatively expects a narrow dust belt composed of 
small grains to be present. While a spectral slope steeper than 
2 is common for debris disks at (sub-)mm wavelengths, such 
a behavior is very unusual in the range of * 100/im among the 
few disks with characterized slope in this wavelength range (e.g.. 



Hillenbrand et al.|2008[|RoccatagUata et al.|2009i ). For the debris 
disks considered there, typical values of A at wavelengths around 
lOOyum are in the range of to 2. Note that our Herschel obser- 
vations are the first to characterize the spectral slope of a large 
number of debris disks at wavelengths around 100 yum. Earlier 
observations in this wavelength range suffer from low sensitiv- 
ity. Thus they were limited to few bright disks or disks with very 
flat spectral slope in the far-ir. The completion of our survey and 
an analysis of the whole data set obtained (Eiroa et al., in prep.) 
will allow to place the three disks presented in a broader context. 
We define an unusually steep SED to be an SED that exhibits 
a Avi.vj > 2 in at least one combination of the PACS bands of 
70yL(m, lOOyum, and 160yum, and for which a Av,_v, = 2 would 
not be possible within the \cr uncertainties. These criteria are 
met by the three debris disks presented in this work (Table|2|l. 



4. SED modeling 

In this section we perform detailed analytical model fitting to the 
observed SED data to explore quantitatively the conclusions on 



the disk properties derived in the above qualitative discussion of 
the steep slope of the observed SEDs. 

4. 1 . Stellar properties 

Table [3] lists some of the main stellar parameters and other rel- 
evant observational properties of HIP 103389, HIP 107350, and 
HIP 114948 (Eiroa et al., in prep.). These parameters were ob- 
tained using the DUNES Virtual Observatory loop] The spec- 



tral types are taken from the Hipparcos catalog (Per ry man et al 
van LeeuwenpOOTJl, and the range of spectral types from 



1997| 

Skiff|( 2010| l. Bolometric luminosities and stellar radii have been 



estimated from the absolute magnitudes and the bolometric 
corrections using the measurements by Flower| ( |199 6); similar 
values are obtained using the bolometric correction procedure 
by Masana et al. (2006). Effective temperatures, gravities and 



metallicities are mean values of spectroscopic and photometric 
estimates from different works in the literature (HIP 103389 and 



HIP 11 4948:, Allsnds Prieto & Lambert, 1999t [Gray et al 2006 



Holmberg et al. 2009[ HIP 107350: Valenti & Fischer 



2005 



Fuhrmann"2008). Radial velocities are taken from 'Mal donado] 
et al. ( 2010 ) and Kharchenko et al. | ^ Wl). Rotational v eloc- 
ities of the stars are taken from Holmberg et al. ( 2007|l and 
Ma rtinez-Arnaiz et am2010). Periods are taken from |Reiners| 
(2006|l and [Messin a et al. (2003). The Li i 6707.8 A line equiva- 
lent widths are taken from Maldonado et al.|(2010| l in the case of 
HIP 107350 and our own estimates using reduced HARPI^data 
for the other two stars; for these latest objects the EW(Li i) have 
been corrected from the Fel 6707.4 A line contamination as in 
Maldonado et al. ( 2010|l. All three stars show log/J^,^ activity 
index levels (|Ciray et al.|:Z006 Martinez- Amaiz et al.|20I0|l in- 



^ http://sdc.cab.inta-csic.es/dunes/searchform.jsp 

"* Based on observations made with the European 



Observatory telescopes 
Archive Facility. 



obtained from the 



Southern 
ESO/ST-ECF Science 
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dicating that they are active stars (e.g., |Maldonado et al.|2010) l. 
The observed ROSAT X-Ray luminosities, Lx/L^oi, derived in 
the context of DUNES (Eiroa et al., in prep.) are also given. 

The masses of the stars in Table [3] are estimated from their 
radii and gravities. We estimate the age of the stars (Table l3]l 
from the log 7?^,^ index, and periods using the relationships as 
given by Mamajek & Hillenbrand ( 2008 ) and from the X-Ray lu- 
minosity following |Garces et al.| ( |2011j ). For all three stars these 
age tracers yield consistent values of less than ~ 500 Myr. We 
note that EW(Li i) values are close to those of the Hyades en- 



For the comparison between modeled and observed measure- 
ments we use a reduced chi-squared 



velope with an age of 600 Myr jRocha-Pinto et al. (2004) give 
an age of 370 Myr for HIP 107350 based on the log^^^ activity 
index, while isochrone-based ages are 2.3 Gyr for HIP 103389, 
with a range of 0.6 Gyr to 3.8 Gyr polmberg et al.|[2009| , 
3.1 Gyr for H IP 107350, with a range of 1.1 Gyr to 4.7 Gyr 
(IValenti & Fischer 2005), and 3.7 Gyr for HIP 114948, with a 



range of 2.1 Gyr to 5.2 Gyr ( [Holmberg et al.||2009| l. Given that 
the stars are located right on the main- sequence, stellar ages are 
difficult to estimate on the basis of isochrones, and that they are 
very sensitive to Teff and metallicity (Holmberg et al. 2009) we 
would favor in these cases the young values estimated from the 
activity tracers and rotational periods. This choice is also com- 
patible with the kinematics of the stars and their ascription to 
the young moving groups, less than 1 Gyr, of HIP 107350 to the 
local association (~ 150 Myr) or the Hercules-Lyra (~ 200 Myr) 
association (Maldonado et al. 2010 Lopez-Santiago et al. 2006 1, 
and HIP 114948 and HIP 103389 to the Castor Moving Group 
(~ 200 Myr) and the young disk population (< 1 Gyr), respec- 
tively. 

HIP 107350 hosts a T dwarf companion (HNPeg B, 
M = 22 Mj) with a projected physical separation of 794 AU 
(«i 43") from the host star ( |Luhman et"aL] |2007|. This is suffi- 
ciently far away, so that the presence of this companion does 
not afffect our present analysis. A very faint source (a; 2 mJy) is 
detected in the PACS 100 yum image close to the predicted po- 
sition of HNPeg B. However, we do not expect HNPeg B to 
be detectable by our observations (» 0.6 mJy at 8 fim). Thus, we 
attribute this detection to a coincidental alignment with a back- 
ground source. 



4.2. Modeling the excess 

In this section, detailed analytical model fftting to the observed 
SED data is performed to explore quantitatively the conclusions 
on the disk properties derived in the above qualitative discussion. 
This is done using the results from Reduction 1, since these re- 
sults give the steep SEDs. The fftting results are expected to dif- 
fer less signiffcantly from standard solutions using the photomet- 
ric data from Reduction 2. For modeling of the three sources, we 
use SAnD which is part of the DUNES modeling toolbox (Lohne 
|et al.|20l"T] Augereau et al., in prep.). A description of this tool 
can be found in the appendix. For the model fftting of the in- 
frared excess SED, the above derived stellar properties are used. 
It is important to note that uncertainties of the stellar parameters 
as well as on the distance of the objects are not considered in 
our modeling of the excess. Thus, the derived uncertainties of 
the disk parameters are only the formal uncertainties from the 
fftting procedure. The uncertainties of the stellar parameters and 
distance would mostly affect the efifective temperature and lu- 
minosity of the star and, thus, the heating of the dust - slightly 
increasing or decreasing its temperature. In particular, they will 
not alTect signiffcantly the slope of the derived grain size distri- 
bution. 



2 
/Vred 






Nm- ^N- Nf,, 



(4) 



where A^ is the number of SED data points used, F,, Fj, and cTi 
are the measured ffux, the modeled ffux, and the uncertainty in a 
data point, A^fiee is the number of free parameters, and A^dot is the 
number of degrees of freedom in our fftting. For this reduced x^ 
a value of 1 .0 is generally desirable in the fftting if A^^of is large, 
while in our case A^dot ranges from 1 to 3 where the reduced x^ 
is desirable to be as small as possible. 

Only selected data points are used. This results in a total of 
8 data points included in the fftting for each source. The SEDs 
of the sources can be seen in Fig. [3] The references for all pho- 
tometric data points available for the three sources can be found 
in Eiroa et al. (in prep.). The selection is done by the following 
criteria: 

- All flux measurements at wavelengths > 10 //m are consid- 
ered, but no upper limits, since they do not give signiffcant 
additional information in the present case. Measurements at 
these wavelengths are also included, if they are photospheric, 
because they exclude signiffcant emission from the disk at 
these wavelengths. 

- Measurements at wavelengths < lOjum are not included, 
since they are at much shorter wavelengths than the shortest 
wavelength at which excess is detected (24yum . . . 32jum). 

- To account for the Spitzer/IRS spectrum in a consistent way, 
an additional photometric point at A- 32 //m is extracted 
from these data following Hillenbrand et al. (2008i and is 
fftted along with the other photometric values. Therefore, 
10 consecutive data points centered at A - 32yum are aver- 
aged. The uncertainty is computed by adding in quadrature 
the standard deviation of these data points and a 5% cali- 
bration uncertainty (Spitzer/IRS instrument handbook v4.0; 
Teplitz et al.|201 1\. This will prove to be sufficient to get the 



whole Spitzer/IRS spectrum properly reproduced by all the 
best-fft models. 

The capabilities of SAnD are used to explore a very broad 
range of parameters. The following model is employed: 

- Radial surface density distribution E (r) oc r^" with inner and 
outer cut-off radii ri,, and rout, 

- Differential grain size distribution dn(a) oc a^'^da with lower 
and upper cut-off size a^i^ and Omax, 

- Two possible grain compositions - pure astronomical sili- 
cate ( |Draine| 2003|l and a 1:1 m ixture of astronomical sili- 
cate and ice ( Lohne et al.|2011 Augereau et al., in prep.) - 
to explore the possibility that water ice might be a signiffcant 
constituent of debris disk dust and might be responsible for 
the peculiar shape of the SEDs. 

Since the number of free parameters in the fftting is close 
to the number of data points used, the fftting is expected to be 
very degenerate. To explore the parameter space in an efficient 
way and to ffnd parameters that can be ffxed, because they have a 
unique best-fft result or they have no signiffcant effect on the fft- 
ting at all, we use different approaches. These approaches are de- 
scribed below. The parameter space explored for each approach 
is listed in Table |4] The fftting results are compiled in Table l5] 
Simulated SEDs from our best-fft models are shown in Fig. [3| 
The parameter space explored and the parameters ffxed in each 
approach are motivated by the results of the previous approaches 
as described in the following. 
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Table 4. Explored parameter space for the fitting based on simulated thermal annealing (for details see Sect. 



4.2 1 



Approach 1 



Approach 2 



Approach 3 



Approach 4 



Parameter 


Range 


# values 


Range 


# values 


Range 


# values 


Range 


# values 


distribution 


^in [AU] 


3-100 


817 


3-100 


817 


3-100 


817 


3-100 


817 


temp 


Rout [AU] 


5-300 


671 


5-300 


671 


5-300 


671 


5-300 


671 


temp 


a 


-2.0-3.5 


56 


fixed 


1 


fixed 


1 


fixed 


1 


linear 


Qmin [yUm] 


0.2 - 20.0 


448 


0.2 - 20.0 


448 


0.2 - 20.0 


448 


0.2-20.0 


448 


logarithmic 


"max [;^rn] 


fixed 


1 


fixed 


1 


fixed 


1 


1.0-1000.0 


674 


logarithmic 


r 


2.0-5.0 


31 


2.0-10.0 


81 


fixed 


1 


fixed 


1 


linear 


A^dust [Msl 


free 




free 




free 




free 




continuous 


Composition 


0%, 50% ice 


2 


0%, 50% ice 


2 


0%, 50% ice 


2 


0%, 50% ice 


2 





4.2.1 . Approach 1 : Exploring the parameter space 

In a first approach, a range of free parameters is considered 
that includes those used to explain most of the known debris 
disks. An upper grain size of 1 mm is chosen, large enough to 
consider any effect of large grains on the fluxes at all wave- 
lengths observed. Seven free parameters are used (nn, rout, «, 
flmin, T, -^dust, ^nd dust chemical composition). This is the largest 
number possible considering a total of 8 SED measurements. 
Although strong degeneracies are expected, this approach is used 
to explore the parameter space without strong initial constraints 
that might bias the results. Parameters that have a unique solu- 
tion or that have no significant effect on the fit can be fixed in 
subsequent approaches. 

A very narrow ring structure is found to be the best-fit for all 
three disks, although the parameters describing the spatial dust 
distribution are not constrained very well. Furthermore, a very 
steep grain size distribution with a value for the exponent y very 
close to the edge of the explored parameter space (5.0) is found 
as best-fit for all three SEDs (Table |5]l. In the following, it is re- 
ferred to this as an underabundance of large grains. Since this 
is an atypical, but not unexpected result (Sect.l3]l, more effort is 
put on the evaluation of the significance of this result in the sub- 
sequent approaches. As the radial extent of the disks has been 
found to be narrow, the exponent of the radial surface density 
distribution is fixed in the further fits to a = 0.0 (constant sur- 
face density), decreasing the explored parameter space by one 
dimension. In a naiTow ring, this will not have any effect on the 
SED. 



4.2.2. Approach 2: A larger range of possible values for y 

Now, the aim is to obtain a best-fit value for y that is included 
in the considered parameter space. A range of possible values 
for y of 2.0 to 10.0 is explored. All other parameters (beside a 
which is now fixed as described above) have the same ranges as 
in Approach 1. This results in a total of 6 free parameters. 

For HIP 107350, the range of explored values of y seems to 
be still too small. This is ignored, since such a very large value 
is not significant (confidence levels: 5.9 . . . 10.0, Table l5]l. For 
the other two objects, the value of y is well within the explored 
range. The confidence levels of the inner and outer radius sug- 
gest that the ring-like shape of the disk is not very significant, 
even for a constant surface density. 



4.2.3. Approach 3: Fixing y to 3.5 

Fixing the value of y to 3.5 and having a fixed upper grain size 
of Qmax = 10 mm, one can force a fit where the parameters of 



the grain size distribution are consistent with an equilibrium col- 
lisional cascade (Dohnanyi 1969). The explored ranges of all 
other parameters are unchanged compared to Approach 2. This 
results in a total of 5 free parameters. 

The resulting x^ is much worse (by a factor of 2.5 to 12.6, 
Table l5]l compared to the results from Approach 2. The changes 
of the radial position of the dust ring are in line with the expecta- 
tions from the much lower abundance of small (warm) particles. 



4.2.4. Approach 4: Fixed value of y 
grain size 



3.5, but free upper 



Another possibility to produce an underabundance of large 
grains in the model is to let the upper grain size be a free param- 
eter Since a very steep size distribution and a small upper grain 
size have comparable results (removing large grains from the 
model), one can fix the value of y to 3.5 to be consistent with an 
equilibrium coUisional cascade. The resulting upper grain size is 
expected to be sufficiently small. This results in a total of 6 free 
parameters. 

This approach gives in general the best;^^^ (Table pb. Upper 
grain sizes of few tens of micron are found. The low radial ex- 
tend of the debris ring is again not very significant, but still the 
best-fit. 



4.2.5. Grid search and inclusion of spatial information for 
HIP 114948 

From our discussion in Sect. |3]we are not able to constrain our 
parameter space by any reasonable physical assumptions. Since 
SAnD is designed to explore a large, high-dimensional parameter 
space using a statistical approach, it was particularly useful to 
find the global best-fit to the data in the above approaches. The 
over-all result of our fitting (a lack of large grains) is clearly vis- 
ible for HIP 103389 and HIP 107350. For HIP 1 14948 the result 
is not that stringent. The value of y = 4.7;^^ 3 from Approach 2 of 
the SAnD fitting is significantly smaller than the values derived 
from the other two disks (although still significantly larger than 
the value of 3.5 expected from a standard equilibrium collisional 
cascade, Dohnanyi 1969). Thus, we perform a grid search of the 
parameter space using GRaTer ( Augereau et al. 1999 ; Lebreton 
et al. 2011^ to confirm our results from the statistical approach 
of SAnD. Furthermore, the constraint that the source is spatially 
unresolved is included in this fitting. This is not expected to 
significantly change the best-fit results (since these results are 
consistent with an unresolved disk) but to further constrain the 
range of values possible within the uncertainties of the derived 
values (since very extended configurations that would clearly be 
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Fig. 3. Observed and modeled SEDs of the three sources. The SpitzerflRS spectrum is plotted in red. The data points considered for 
the fitting are shown in dark blue. The PACS data considered are from Reduction 1, while the results from Reduction 2 are plotted 
in light blue for comparison. Top: Final models from the SAnD approaches. The modeled SEDs (solid line) for the three disks are 
computed from the results of Approach 4 using pure astronomical silicate (Sect. |5]). The dashed line in each panel represents the 
model stellar photosphere. Bottom: Comparison between the best-fit results from the different approaches. The modeled SEDs are 
computed from the best-fit result of each approach using pure astronomical silicate. The results from Approach 1 are not included, 
since the best-fit results are found not to be in the explored range of parameters in this approach. 



resolved are also included in the confidence levels from pure 
SED fitting). The results from fitting with and without the spa- 
tial constraints are compared. 

Since the general result from the SAnD fitting (steep grain 
size distribution, large lower grain size) is independent from 
the approach used, we decide to use the standard approach of 
GRaTer This is a power-law approach (dn (a) oc a^'^da) with up- 
per and lower limit for the differential grain size distribution and 
a two-power-law approach for the radial surface density distri- 
bution; 



nR) ^ j2{iR/Ror'-"'' + (R/Rof''y 



(5) 



where a and «!„ are slopes of the radial surface density distribu- 
tion outside and inside a peak position. We fix ain = 10 to get a 
sharp inner edge close to Rq comparable to the single power-law 
approach used for the SAnD fitting. 

The spatial constraints are included in our fitting using a ra- 
dial brightness profile derived by radial averaging of the sources 
image in our lOOjum data. After a Gaussian is fitted to the image 
to derive the source center, the average flux and distance from 
the center is derived for all pixels in one pixel wide radial bins. 
The uncertainty on the flux in each bin is derived using the stan- 
dard deviation of the pixel values in each radial bin. Images are 



simulated from the models, convolved with the PACS PSF mea- 
sured on the bright standard star a Boo and profiles are extracted. 
These profiles are compared to the one derived from the data. 
The same number of SED data points as before is considered. In 
addition, a total of 10 profile points are used. An over-all red uced 
X^ including SED and profile data is computed following Eq, 



A.l 



(using only one profile). The explored parameter space and best- 
fit values for the fit with and without the profile included are 
given in Table 6] These results should be compared to the results 
from Approach 2 of the SAnD fitting, since these two approaches 
are the most similar They are found to be consistent. Maps of the 
distribution of the reduced x^ in the parameter space for the fits 
including and not including the profile have been produced as 
cuts through the parameter space in three planes. The remaining 
parameters are fixed to their best-fit values, respectively. These 
maps are displayed in Fig. HI We find that the best-fit values do 
not change significantly when including the spatial information. 
This confirms that the results from the pure SED fitting are fully 
consistent with a spatially unresolved disk. Fig. [4] illustrates a 
number of results from the mapping of the parameter space: 

- There is a correlation between the lower dust grain size 
Omin and the exponent -y of the size distribution. However, 
from our SED fitting one can rule out that any combination 
of significantly smaller grains and smaller exponent of the 
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Table 6. Fitting result for HIP 1 14948 using the grid search method with GRaTer 



Parameter 



Range 



Spacing 



# values 



Fit w/o spatial constraints 



Fit w/ spatial constraints 



Ra [AU] 




1.0-250.0 




logarithmic 


60 


13.7 


a 




0.0-10.0 




linear 


11 


8.0 


flniin [yUm] 




0.05-43.5 




logarithmic 


55 


14.1 


r 




2.7-5.9 




linear 


17 


4.7 


V,ce/(Vs. + Vice) 




0.0-0.9 




linear 


10 


0.6 


Md„st [Me] 


1.0 X 10-*- 4.0 X 


10-^ 


logarithmic 


100 


2.36 X 10-5 


reduced x^ 




- 




- 


- 


0.19 



13.7 
9.0 
14.1 
4.7 
0.6 
2.20 X 10"' 
0.12 




100.0- 





nSO.O 



40.0 



10.0 



5 4 3 5 4 3 0.1 1.0, , 10.0 




100.0 



10.0 





30.0>E 



20.0 o 



10.0 
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)' }' «min[lim] 



0.0 



Fig. 4.x^ maps for HIP 1 14948. The maps are created as two di- 
mensional cuts through the;i^^ distribution in the parameter space 
searched with GRaTer without (top) and with {bottom) the spa- 
tial information included. The parameters not plotted are fixed 
to their best-fit value, respectively (Table [6]). The correlations 
between the three most relevant parameters (y, flmin, and Rq) are 
shown. The white cross denotes the position of the best-fit (low- 
est x^)- The dashed vertical line denotes the value of y = 3.5 
expected from a standard equihbrium collisional cascade. 



size distribution results in reasonable fits to the data (which 
would be closer to the expectations from a collisional cas- 
cade and grains close to the blow-out size of the system of 
a < l.Qfim). 

- There is nearly no coiTelation between the radial position of 
the dust belt and the exponent of the grain size distribution 
in both the case without and with the spatial information in- 
cluded. 

- The degeneracy between the lower dust grain size and the 
radial position of the dust belt known from SED fitting is 
clearly visible. It is also obvious that the inclusion of the 
spatial information in part breaks this degeneracy since it ex- 
cludes a large (> 70 AU) radial distance of the belt from the 
central star 

- In particular, the results of a^i^ and y seem to be less con- 
strained when including the spatial information. This is in 



part because the spatial information gives no constraints on 
the parameters of the grain size distribution (beside the partly 
broken degeneracy of the lower dust grain size and the radial 
position of the dust which results in an indirect constraint on 
flmin)- Thus, including the radial profile we include a num- 
ber of measurements that are not sensitive to changes in y, 
which results in a less significant increase of the reduced x^ 
with changes of y in this case. In addition, the data points in 
the profile can not be considered to be independent from each 
othe r. Alt ho ugh t his is mitigated in part by the weighting in 



Eqs.|A.l|to 



A.3 



the combination of ;t^^ derived from SED 
and profile is not fully consistent, in particular in the case of 
unresolved sources. To measure the quantity and uncertainty 
of a model parameter through model fitting using reduced 
X^, one should in general only include measurements in the 
fitting process that are sensitive to the parameter (as far as 
one can tell in advance) and independent from each other. 
Thus, the uncertainties on y derived from pure SED fitting 
should be taken more reliable here. 



5. Results 

The general results and conclusions that can be drawn from the 
fitting approaches are summarized in the following. 

- Values of y = 4.7 to 10.0 or a small upper grain size are 
found rather than the expected value of y x; 3.5 and a large 
«max- Models with y = 3.5 do not reproduce the observed 
SEDs in a reasonable way. This result is not strongly affected 
by any modeling degeneracies found (that are included in the 
estimate of the uncertainties). 

- Evidence is found for a large lower grain size of the 
dust compared to the expected blow-out size of the sys- 
tems (0.5 yum to 0.9 yum). However, the lower grain size 
depends very much on the approach used. Earlier studies 
of debris disk SEDs had similar results (e.g., HD 107146; 
Roccatagliata et al.]|2009 1 that have been at least mitigated 



and partially attributed to modeling degeneracies in pure 
SED fitting after including resolved data in the fitting pro- 
cess ( jErtel et al.||2011| ). Thus, this result has to be treated 
with caution. 

- The dust in all three disks must be located at a distance from 
the star of no more than a few tens of AU as can be con- 
cluded from the modeling results and the fact that the disks 
are spatially unresolved in the PACS images. Stronger con- 
straints on this value are not possible due to degeneracies in 
the modeling. 

- The disks appear to be narrow rings. However, the uncertain- 
ties - in particular of the outer disk radius - are very large 
and also very broad disks can result in fits on the SED that 
fall within the derived 3 cr confidence levels. This can only 



S. Ertel et al.: A peculiar class of debris disks from Herschel/DUNES 



partly be excluded by the additional constraint that the disks 
are spatially unresolved. 

It is important to note that these results depend in part on the 
model used and on the parameter space explored. For example, 
the composition and shape of the dust grains have not been ex- 
plored in detail in our modeling, but might allow one to repro- 
duce the SEDs with models that are more in Une with a standard 
equilibrium collisional cascade. 

Furthermore, the results from Reduction 2 have not been 
considered in the fitting. Considering these photometric mea- 
surements and uncertainties instead of those from Reduction 1 
would lower the significance of the underabundance of large 
grains. As a result, values are expected to be found for the model 
parameters that are more in line with those expected from an 
equilibrium collisional cascade. Without a deeper understanding 
of the data obtained with Herschel, it is not possible to make a 
final statement about the significance of our modeling results. 

Each approach presented results in a different set of best-fit 
parameters and a different ;t'^, while the results can be interpreted 
in a consistent way. To find from the fitting results a final model 
of each disk, the following selection is applied: 

- The results from Approach 1 have to be ruled out, because 
the best-fit parameters from this approach are close to the 
edge of the explored parameter space. 

- In general, the SEDs cannot be fitted by a grain size dis- 
tribution that includes a significant amount of large grains 
(Approach 3). 

- The results from Approach 2 (free y, a^ax = 1.0 mm) and 
Approach 4 (free amax, T = 3.5) are consistent, while the re- 
sults from Approach 4 and pure silicate represent in general 
the best fit to the data. Thus, these results are considered as 
the best-fit models from the fitting. 

From these final models, one finds dominating dust temperatures 
(the temperature of the smallest grains that dominate the emis- 
sion) of 63 K, 50 K, and 73 K for HIP 103389, HIP 107350, and 
HIP 1 14948. 



5. 1 . The origin of tlie dust 

The existence of micron-sized dust grains around main sequence 
stars is usually explained by collisions of larger bodies such as 
planetesimals (e.g., Krivov|2010 i. In such a scenario one would 
expect a significant amount of larger grains (several tens of mi- 
cron to millimeter sized particles) to be present that are pro- 
duced by the collisional cascade also producing the smallest 
grains ( Dohnanyifl 969) . These grains significantly contribute to 
the SED in the observed wavelength regime. Assuming that the 
steep slopes of the SEDs are real, we find from our modeling 
strong evidence for an underabundance of these larger grains. 
To our knowledge, these are the first debris disks discovered, 
that potentially exhibit such a peculiar shape of the SED in the 
wavelength regime of 70yum to 160//m. All this would sug- 
gest that these objects are exceptions from the common under- 
standing of dust creation in debris disks. On the other hand, six 
more disks with similar shape of the SED have recently been 
identified in our DUNES survey, which might suggest that this 
phenomenon is common among low-mass (i.e., low luminosity. 
Tables Is] and ^ debris diskaj The fact that they are very faint 
(Table p| might imply that we are faced with a new class of de- 
bris disks that were not revealed earlier due to limited sensitivity 



in the relevant wavelength range. In the following, we discuss 
scenarios that might be capable to explain the presence of a disk 
as it has been modeled for the three systems. 

Scenario 1: Significant deviation from the conditions re- 
quired for a standard equiUbrium colUsional cascade 

For the standard equilibrium collisional cascade, a number of as- 
sumptions are made that are not necessarily valid in debris disks 
due to the effects of radiation pressure and Poynting-Robertson 
(PR) drag (Sect. [Til. For the massive debris disks known so far, 
PR drag is expected to be negligible (Wyatt 2005). Following 
Backman & Paresc^ ( |1993| , one can estimate the ratio between 
collisional time scale and PR time scale for the three disks 
to 10"^ . . . 10^^. However, the disk models used (best-fit from 
Approach 4, pure silicate) represent very narrow rings, which 
is not a significant fitting result. Assuming a ring with a width 
of 10 AU starting at the inner radius found from the fitting, one 
finds fcou/fpR = 3 X 10"^ ... 5 X 10"^. It is not possible to put 
strong constraints on the dust dynamics due to the crude esti- 
mate of the time scales. Thus, transport mechanisms may play a 
role in the dust dynamics of the small grains in these disks. The 
eff'ects of these processes have been modeled, e.g., by Krivov 



eraL] ( |2006l l; |Thebault & Augereau| ( |2007l ); |Wyatt eraL] ( |2011| l 
In particular, radiation pressure leads to a wavy structure of the 
grain size distribution increasing the slope of the grain size dis- 
tribution in certain ranges. PR drag results in a depletion of small 
grains, increasing the dominating grain size, but mitigating the 
wavy structure of the size distribution. 

Scenario 2: Diff'erent chemical composition or physical shape 
of the dust grains than the assumed one 

One might imagine a dust composition emitting significantly 
more efficiently than astronomical silicate or ice in the wave- 
length range of lOfim to 100/im compared to longer wave- 
lengths. Therefore, one would need grains with a break in 
2abs (^) (Sect. l3]l at shorter wavelengths (in particular for large 
grains). Also, if the large grains were significantly colder than 
expected, their emission would be reduced. This would both al- 
low a significantly larger amount of large grains to be present in 
the systems than modeled. Furthermore, Mie calculations might 
result in an incorrect approximation of Qabs particularly for large 
grains (tens of microns or larger), e.g., due to incorrect assump- 
tions on the physical shape of the grains. 

Voshchinnikov et aT] ( |2006[ ) simulated the effect of poros- 



ity on the absorption efficiency of dust grains. They found that 
the temperature of dust grains at a given distance from a star 
decreases significantly with increasing porosity. If one now has 
small, compact grains and larger, porous grains, the difference 
in temperatures of the smaller and larger grains can be increased 
significantly. This would be the case, if the large grains produc- 
ing the smaller ones through collisions were porous, composed 
of smaller, compact units in the order of the lower grain size 
derived from the modeling of the three systems (i.e., 5//m to 
lOyum). Such a scenario would result in an overabundance of 
grains of this size compared to the larger, weaker grains. Such 
an overabundance would shift the lower cut-off of our adopted 
power-law distribution toward this grain size, although smaller 
grains are present in the system. 

Scenario 3: A shepherding planet 

Planets can trap dust particles into mean-motion resonances 
(e.g., Wyatt 2006 1. This results in a barrier against particles mov- 



^ For an overview of the DUNES results we refer to Eiroa et al., in 
prep. 



ing inward due to PR drag. This barrier is less efficient for very 
small grains for which PR drag is very strong ( Reidemeister et al. 
201 1 1. On the other hand, very large grains and planetesimals are 
not significantly affected by PR drag. 
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Table 7. Simulated observational properties of the disks derived from our final models 



Source Ldi,st/Z,« 



HIP 103389 1.5 X 10^' 
HIP 107350 0.6 X 10^^ 
HIP 114948 2.5 X 10^"* 



(■rdusl/^*)o.6;mi 

face-on edge-on 



1.9 X lO^-* 
7.4 X 10^'^ 
3.4 X 10^" 



(-^du.st/-P'*)o.8^m 

face-on edge-on 



face-on edge-on 



face-on edge-on 



1.1 X 10" 
4.9 X 10^' 
2.1 X 10" 



1.1 X 10" 
4.6 X 10^' 
2.0x10" 



1.7 X 10-^ 
8.5 X 10" 
3.0x10-^ 



1.3x10" 

5.5 X 10-^ 

2.6 X 10" 



1.8 X lO^** 
1.3 X 10" 
3.3 X 10-* 



1.8 X 10" 

6.9 X 10-' 
3.7 X 10" 



1.5 X 10-* 
3.9 X 10" 

2.6 X 10^-* 



The dust seen in these disks might be produced in a faint, 
transport dominated (Krivov 2010[ l debris disk further away 
from the star, too cold and too faint to be detected in the available 
data. The small particles (few tens of micron and smaller) would 
then be dragged inward by PR drag. A possible planet present 
further inward in the system should trap the particles into reso- 
nance. This would result in an accumulation of particles with a 
very distinct range of sizes, which would explain both the lack 
of large grains (at least in an abundance and at a position where 
they would be detectable) and the lack of small grains. Such a 
"dust trail" has been observed, e.g., to be associated with our 
Earth and has been modeled to predominantly consist of grains 
of X llfim and larger (Dermott et al. 1994). It will have to be 
evaluated, whether an outer planetesimal belt capable to produce 
enough dust through collisions to replenish the dust in the disk 
can be faint enough not to be detected by the Herschel obser- 
vations, and under which conditions the planet can trap enough 
dust into resonance to produce the disk observed. This has to be 
done through detailed dynamical and collisional modeling (e.g., 
Reidemeister et al.|2011[ l, which is not within the scope of the 
present work. 

5.2. Observational perspectives 

Further observations with present and near future instruments 
can help to reasonably increase our understanding of these debris 
disks. In the following, we discuss a number of observations that 
might result in very valuable data. 

- Since the most extended models of the three disks predict 
radial extents of only » 2'.'2, these objects are at the edge 
of what is resolvable with coronagraphy . Furthermore, the 
disks are very faint in scattered light (Table ITli. It is impor- 
tant to note that in the case of edge-on orientation most of the 
flux comes from forward scattering and will then be concen- 
trated close to the star (projected physical separation) with 
no contribution to the signal in coronagraphic observations. 
A contrast ratio of < 10"'' is not accessible to present instru- 
ments. 

- Successful optical/near-infrared imaging of planetary com- 
panions would give strong evidence that Scenario 3 (shep- 
herding planet) is responsible for the peculiarities of the disk. 
Determining the position of this planet would also help to 
further constrain the position of the dust in this scenario. The 
youth of the stars means good chances to directly image gi- 
ant planets at separations of > 1" from the star with present 
methods (e.g., Marois et al. 2010). 

- The predicted extents of the disks are in a range easily resolv- 
able with ALMA. However, the disks are already very faint 
at PACS wavelengths and the expected surface brightness 
is decreasing particularly steep toward wavelengths acces- 
sible with ALMA. It is not clear, whether ALMA observa- 
tions are sensitive enough to detect these disks. On the other 



hand, at (sub-)mm wavelengths even upper limits in the sen- 
sitivity range reachable are expected to provide useful fur- 
ther constraints on the shape of the SEDs. From the results 
of our fitting Approach 4, predicted fluxes of star and disk 
at 350yL(m (the shortest wave band accessible with ALMA) 
are 0.55 mJy and 0.11 mJy for HIP 103389, 0.53 mJy and 
0.05 mJy for HIP 107350, and 0.60 mJy and 0.20 mJy for 
HIP 114948, respectively. 
- Further photometry and spectroscopy in thermal reemission 
would help to constrain the radial distribution of the dust 
and the size distribution of the small grains (wavelengths 
between 40 //m and 70;um) and the steep slope and the 
shape of the SEDs at longer wavelengths {A > 70//m). The 
Stratospheric Observatory for Infrared Astronomy (SOFIA) 
is most promising to provide observational capabilities in the 
relevant wavelength regime in the near future. 

6. Conclusions 

The first data at wavelengths > 70;um for the three debris disks 
discussed in this paper have been presented and modeled. All 
three sources potentially exhibit an unusually steep decrease of 
the SED in the wavelength range between 10 pm and 160 ;um. In 
a general discussion, it has been shown that this peculiar shape 
of the SED is an indicator for a deviation from the case of a stan- 
dard equilibrium collisional cascade assuming standard grain 
composition (astronomical silicate or ice) and shape (compact, 
spherical grains). However, it has also been shown that the un- 
derstanding of faint source photometry obtained with Herschel 
is still incomplete and that the results presented here depend 
very much on the actual version of the data reduction pipeline. 
Provided that the steep decrease of the SEDs is real, modeling 
implies that the thermal emission from these disks is dominated 
by a very distinct grain size regime of several micron to few 
tens of micron. The disks have been modeled as narrow rings 
with a significant underabundance of large grains. A number of 
possible explanations for such an unusually steep shape of the 
SEDs have been discussed. Six more candidates for this new 
class of debris disks have been identified so far from the on- 
going DUNES survey. This is the first published discovery of 
debris disks that potentially exhibit such a peculiar shape of the 
SED in this wavelength regime. 
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Appendix A: Description of SAnD 

When performing analytical model fitting of debris disks, one 
is faced with a number of challenges. On the one hand, re- 
cent results (e.g., [Backman et al.|2009 Ertel et al.|2011 1 reveal 
more and more the complexity of the known debris disks and 
the inability to reproduce the available data with simple mod- 
els. Furthermore, standard assumptions (e.g., lower grain size 
consistent with the blow-out size of the system, grain size dis- 
tribution exponent of -3.5) have been proven to be inadequate 
(e.g., Ertel et aL]|2011 this work), and unexpected results like 
an outward increasing surface density distribution (Lohne et al.| 



2011 1 have been found. The increasing number and quaUty of 



the data available (in particular, spatially resolved images from 
Herschel and ALMA) allow one to break modeling degenera- 
cies through simultaneous multi-wavelength modeling and to in- 
crease the complexity of the models. On the other hand, finding 
a best-fit model in a complex parameter space results in signif- 
icant challenges on the fitting method and the simulation of the 
observed data in order to reduce the computational effort neces- 
sary. 

To rise to these challenges we, developed the tool SAnD. It is 
able to fit SED data and radial profiles in thermal reemission si- 
multaneously using a simulated annealing approach ( [Press et ST| 
ri992) on a grid of possible values in the parameter space. It uses 
an analytical approach for the radial density distribution of a ro- 
tationally symmetrical disk as well as for the grain size distri- 
bution. The SED is computed in an analogous way to Wolf & 
[Hillenbrand' (2005). Resolved images are computed with suffi- 
cient resolution and convolved with the telescope PSF. Radial 
profiles are then extracted along the major and/or minor axis. 
Each image (i.e., each pair of major and minor axis profiles or 
single major axis profile, if no minor axis profile is extracted, 
e.g., in azimuthally symmetrical disk images) is scaled by a fac- 
tor xi to minimize the;^^ ^^^ (Eq. A. 3 1. This way, only the shape 
of a radial profile is fitted. This is done because flux calibration 
uncertainties have to be considered only once per wavelength, 
but would be included for each data point in the radial profiles 
if absolute profiles were considered. Furthermore, deficits of our 
models to reproduce the absolute flux at one wavelength at which 
also radial profiles are fitted (e.g., due to uncertainties in the op- 
tical properties of the dust) would be considered in each profile 
data point and the SED data point which would result in an over- 
weighting of the data at this wavelength. Simulated and observed 
multi-wavelength data are then compared using the following re- 
duced ;i'^: 






with xlzD = Yj 

2 N:j 

and xl,rof = Z Z 



Fk 



o-k 



-"'i'J i,i,k 'J i,j,k 



(A.l) 
(A.2) 
(A.3) 



where / = is for the SED and / > 1 for the profiles, y = 1 for 
major axis and j -2 for minor axis. In equations A.l to A.3 the 
notations used are: 
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- M - number of images 
-No- number of SED data points 

- A^,-,i - number of profile data points along major axis 

- M,2 - number of profile data points along minor axis 

M 2 

- A^(ot - total number of data points, A^o + 2 2 Nj j 

'■=1 ;=i 

- A^free - number of free parameters 

- Fj., Fi; - modeled and observed SED, respectively 

- S ij^k, S ij^k - modeled and observed surface brightness pro- 
files 

- (Tk, cTijj. - observed uncertainties 

- Xi - surface brightness profile scaling factor 

- Wi - weight. 

For the present work, only one profile is used and wq and 
Wi are set to 1.0, which means that the profile is given the same 
weight as the SED. 

In order to allow one to find the best-fit model within a given 
range of parameters of a high dimensional parameter space in 
reasonable time, one has to compute the value of ;^^ for a given 
set of parameters very fast. Creating the SED and three profiles 
at different wavelengths for one particular model and comparing 
them with the observations usually takes ~ 0.1 sec on an Intel 
Xeon E5410 CPU (2.3 GHz). The exact time depends mostly on 
the particular set of parameters and the number of wavelengths, 
at which radial profiles are provided. It hardly depends on the 
number of data points in the SED and the profiles. 

The simulated annealing approach enables us to handle a 
very large and high dimensional parameter space without sam- 
pling the whole range of parameters. Therefore, a random walk 
on the grid of possible parameters is started at an arbitrary posi- 
tion. The probability p to go a certain step depends on x\ of the 
actual and Xq of the the previous step. It is 1 for xl ^ X^ (fit be- 
comes better) and follows a Boltzmann distribution for Xq < x\ 
(fit becomes worse): 



p = exp 



2 2 



(A.4) 



In this equation, the quantity T is a parameter comparable to the 
Temperature of an "annealing" physical system and controls the 
ability to reach areas in the parameter space that give worse fits. 
It is chosen at the beginning of each run to be a large value based 
on experiences from previous runs and is lowered with each suc- 
cessful step of the random walk. The width of each step is chosen 
from a probability distribution that prefers short steps, while the 
maximum step width is 1 /2 of the whole range of values for each 
parameter. The probability to go long steps is lowered with each 
successful step. With this approach the code is able to reach each 
position in the parameter space, but prefers regions with low;^'^. 
The run stops, when a position in the parameter space with the 
following properties is reached: 

- With the actual value of T the code was not able to leave the 
position after computing 500 further models. 

- The actual x^ is lower than or equal to the best reached x^ 
during the whole run. 

- The actual ;i'^ is lower than or equal to a maximum ;^^^ based 
on experiences from previous runs (input parameter). 

A set of model parameters that satisfies all of these conditions is 
considered as a likely global fit. To increase confidence in this 
result, one can start several runs at different regions of the pa- 
rameter space (4 in the recent case for the final fitting). Only if 



all runs reach the same region in the parameter space, we con- 
sider the result to be really a global best-fit. If only the first and 
not more than one other of these stopping criteria is satisfied, the 
value of T is increased and the code goes on trying to leave the 
actual position in the parameter space. 

Error estimates are done in the context of the simulated an- 
nealing approach by starting a new random walk at the best-fit 
position with a fixed value of T (usually 10% of the best-fit x^ 
in our case). The code counts how often each value of a param- 
eter is reached. This way, the code samples the projection of the 
probability distribution of the random walk (Eq. |A.4[ ) on the axis 
of each parameter With the knowledge of the value of T one can 
then compute the probability distribution 



p' =exp[-{x^-xL,)] 



(A.5) 



and corresponding levels of confidence. 

As a feature of the simulated annealing approach the time 
it takes to find a best-fit (and uncertainties) does not vary very 
much with the number of free parameters or the total size of the 
parameter space to be searched, but depends on the shape of the 
X^ distribution in the parameter space. Weak variations of the x^ 
over a broad range of values of the parameters as well as very 
degenerate fitting problems require more time, i.e. more mod- 
els to be computed. For the recent work (SED fitting only) we 
find the best-fit in a parameter space composed of ~ 10'" to 10" 
knots typically in less than 24 hours including the estimates of 
the uncertainties. Therefore the code computes about 100.000 to 
1.000.000 models and accepts (goes) about 10.000 and 100.000 
steps for the fitting and computes ~ 1.000.000 models and ac- 
cepts ~300.000 steps for the error estimates. 
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